R Notebook

##Loading necessary packages and data

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
Warning: package ‘ggplot2’ was built under R version 4.2.1Warning: package ‘dplyr’ was built under R version 4.2.1── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(janitor)

Warning: package ‘janitor’ was built under R version 4.2.1
Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

    chisq.test, fisher.test

library(lubridate)

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

hb_agesex <- read_csv("../phs_scotland_group_project/raw_data/covid_raw_data/hospital_admissions_hb_agesex_20220302.csv") %>% 
  clean_names

Rows: 43516 Columns: 12── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (8): HB, HBQF, AgeGroup, AgeGroupQF, Sex, SexQF, AdmissionType, AdmissionTypeQF
dbl (4): WeekEnding, NumberAdmissions, Average20182019, PercentVariation
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

board_names <- read_csv("raw_data/covid_raw_data/health_board_names.csv") %>% 
  clean_names

Rows: 18 Columns: 5── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (3): HB, HBName, Country
dbl (2): HBDateEnacted, HBDateArchived
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#data exploration

hb_agesex

board_names

hb_agesex <- left_join(hb_agesex, board_names, by = "hb")

hb_agesex <- hb_agesex %>% 
  select(week_ending, hb, hbqf, age_group, age_group_qf, sex, sex_qf, admission_type, admission_type_qf, number_admissions, average20182019, percent_variation, hb_name.x) 

Error in `select()`:
! Can't subset columns that don't exist.
✖ Column `hb_name.x` doesn't exist.
Backtrace:
  1. hb_agesex %>% ...
  3. dplyr:::select.data.frame(...)
  6. tidyselect::eval_select(expr(c(...)), .data)
  7. tidyselect:::eval_select_impl(...)
 16. tidyselect:::vars_select_eval(...)
     ...
 19. tidyselect:::reduce_sels(node, data_mask, context_mask, init = init)
 20. tidyselect:::walk_data_tree(new, data_mask, context_mask)
 21. tidyselect:::as_indices_sel_impl(...)
 22. tidyselect:::as_indices_impl(x, vars, call = call, strict = strict)
 23. tidyselect:::chr_as_locations(x, vars, call = call)

##the data is relatively organised, let’s start creating some plots

##demographic investigation

hb_simd <- read_csv("raw_data/covid_raw_data/hospital_admissions_hb_simd_20220302.csv") %>% 
  clean_names()

Rows: 21138 Columns: 9── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (4): HB, HBQF, AdmissionType, AdmissionTypeQF
dbl (5): WeekEnding, SIMDQuintile, NumberAdmissions, Average20182019, PercentVariation
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

hb_simd

hb_simd <- left_join(hb_simd, board_names, by = "hb")

hb_simd <- hb_simd %>% 
  select(week_ending, hb, hbqf, simd_quintile, admission_type, admission_type_qf, number_admissions, average20182019, percent_variation, hb_name) 
 hb_simd

hb_simd <- hb_simd %>% 
  mutate(hb_name = if_else(
    is.na(hb_name),
    "Scotland",
    hb_name
  ))

hb_simd <- hb_simd %>%
  mutate(week_ending = ymd(week_ending))

hb_simd <- hb_simd %>% 
  mutate(month = month(week_ending, label = TRUE),
         year = year(week_ending))
hb_simd

hb_simd %>% 
  distinct(simd_quintile)

hb_simd %>% 
  filter(year == "2020", 
         hb_name == "Scotland", 
         admission_type == "Emergency") %>% 
ggplot(aes(x = simd_quintile, y = number_admissions))+
  geom_col()

hb_simd %>% 
   filter(admission_type == "Emergency")%>% 
  group_by(hb_name, week_ending, simd_quintile) %>% 
  summarise(mean_admissions = mean(number_admissions),
            mean_20182019_admissions = mean(average20182019)) %>% 
   
  ggplot() +
    geom_line(aes(x = week_ending,
                  y = mean_admissions), group = 1, colour = "red") +
  geom_line(aes(x = week_ending,
                y = mean_20182019_admissions), group = 1, colour = "blue") +
  scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size =7)) +
   geom_vline(xintercept = as.numeric(as.Date("2020-01-01")), linetype=4)+
  geom_vline(xintercept = as.numeric(as.Date("2021-01-01")), linetype=4)+
  geom_vline(xintercept = as.numeric(as.Date("2022-01-01")), linetype=4)+
  facet_wrap(~hb_name, scales = "free_y") + 
  labs(x = "Date",
       y = "Mean admissions per week",
       colour = "Covid vs Pre-Covid")

`summarise()` has grouped output by 'hb_name', 'week_ending'. You can override using the `.groups` argument.

hb_simd %>% 
   filter(admission_type == "Emergency", 
          simd_quintile == "1",
          year == "2022")%>% 
  group_by(hb_name, week_ending, simd_quintile, year) %>% 
  summarise(mean_admissions = mean(number_admissions),
            mean_20182019_admissions = mean(average20182019)) %>% 
  
  ggplot(aes(x = week_ending, y = mean_admissions))+
  geom_col()

`summarise()` has grouped output by 'hb_name', 'week_ending', 'simd_quintile'. You can override using the `.groups` argument.

hb_simd

hb_simd %>% 
   filter(admission_type == "Emergency", 
          simd_quintile == "1")%>% 
  group_by(hb_name, week_ending, simd_quintile) %>% 
  summarise(mean_admissions = mean(number_admissions),
            mean_20182019_admissions = mean(average20182019)) %>% 
  
  ggplot(aes(x = week_ending, y = mean_20182019_admissions))+
  geom_col()

`summarise()` has grouped output by 'hb_name', 'week_ending'. You can override using the `.groups` argument.

#simd 1 2020
hb_simd %>% 
   filter(admission_type == "Emergency", 
          simd_quintile == "1", 
          year == "2020", 
          hb_name == "NHS Lothian")%>% 
  group_by(hb_name, week_ending, simd_quintile) %>% 
  summarise(mean_admissions = mean(number_admissions),
            mean_20182019_admissions = mean(average20182019)) %>% 
  
  ggplot(aes(x = week_ending, y = mean_admissions, colour = simd_quintile))+
  geom_line()+
   scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") 

`summarise()` has grouped output by 'hb_name', 'week_ending'. You can override using the `.groups` argument.

hb_simd$simd_quintile <- factor(hb_simd$simd_quintile, levels = c("1", "2", "3", "4", "5"))
hb_simd

#simd 5
hb_simd %>% 
   filter(admission_type == "Emergency")%>% 
  group_by(hb_name, week_ending, simd_quintile) %>% 
  summarise(mean_admissions = mean(number_admissions),
            mean_20182019_admissions = mean(average20182019)) %>% 
  
  ggplot(aes(x = week_ending, y = mean_20182019_admissions, colour = simd_quintile))+
  geom_line()+
   scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") 

`summarise()` has grouped output by 'hb_name', 'week_ending'. You can override using the `.groups` argument.

hb_simd %>% 
   filter(admission_type == "Emergency", 
          hb_name == "Scotland") %>% 
  group_by(hb_name, week_ending, simd_quintile) %>% 
  summarise(mean_admissions = mean(number_admissions)) %>% 
  ggplot()+
  geom_line(aes(x = week_ending, y = mean_admissions, colour = simd_quintile))+
   scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") 

`summarise()` has grouped output by 'hb_name', 'week_ending'. You can override using the `.groups` argument.

hb_simd %>% 
 filter(admission_type == "Emergency") %>% 

  
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(lng = -4.975, 
                   lat = 55.445, 
                   color = "blue",
                   popup="Ayrshire and Arran",
                   radius = sqrt(98202/ 10), weight = 1) %>% 
  addCircleMarkers(lng = -2.83333000, 
                   lat = 55.58333000, 
                   color = "blue",
                   popup="Borders",
                   radius = sqrt(20656/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.857784, 
                   lat = 54.988285, 
                   color = "blue",
                   popup="Dumfries and Galloway",
                   radius = sqrt(33981/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.78535, 
                   lat = 56.0021, 
                   color = "blue",
                   popup="Forth Valley",
                   radius = sqrt(64750/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -2.988, 
                   lat = 57.228, 
                   color = "blue",
                   popup="Grampian",
                   radius = sqrt(108032/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -4.71, 
                   lat = 57.12, 
                   color = "blue",
                   popup="Highland",
                   radius = sqrt(60411/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.083999664, 
                   lat = 55.905496378, 
                   color = "blue",
                   popup="Lothian",
                   radius = sqrt(185589 / 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.0, 
                   lat = 59.0, 
                   color = "blue",
                   popup="Orkney",
                   radius = sqrt(3375/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -1.2689, 
                   lat = 60.3038, 
                   color = "blue",
                   popup="Shetland",
                   radius = sqrt(3248/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -7.02, 
                   lat =  57.76, 
                   color = "blue",
                   popup="Western Isles",
                   radius = sqrt(5398/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.1999992, 
                   lat =    56.249999, 
                   color = "blue",
                   popup="Fife",
                   radius = sqrt(84101/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.7333304, 
                   lat = 56.6999972, 
                   color = "blue",
                   popup="Tayside",
                   radius = sqrt(104225/ 10), weight = 1) %>% 
   addCircleMarkers(lng = -4.4057, 
                   lat = 55.90137, 
                   color = "blue",
                   popup="Greater Glasgow and Clyde",
                   radius = sqrt(273189 / 10), weight = 1) %>% 
   addCircleMarkers(lng = -3.83333, 
                   lat = 55.583331, 
                   color = "blue",
                   popup="Lanarkshire",
                   radius = sqrt(143644 / 10), weight = 1) 

NA

hb_simd$simd_quintile <- factor(hb_simd$simd_quintile, levels = c("1", "2", "3", "4", "5"))
hb_simd

library(plotly)

Warning: package ‘plotly’ was built under R version 4.2.1Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(plotly)
hb_simd$simd_quintile <- factor(hb_simd$simd_quintile, levels = c("1", "2", "3", "4", "5"))
hb_simd

#mean admissions per SIMD from 2020 to 2022
p2 <- hb_simd %>% 
   filter(admission_type == "Emergency",
          # hb_name reactive
          hb_name == "Scotland")%>% 
  group_by(hb_name, week_ending, simd_quintile, year) %>% 
  summarise(mean_admissions = mean(number_admissions),
            mean_20182019_admissions = mean(average20182019)) %>% 
  
  ggplot()+
  geom_point(aes(x = week_ending, 
                  y = mean_admissions, colour = simd_quintile, 
                 text = paste0("Date: ", week_ending, 
                               "<br>",
                               "Mean admissions: ", 
                              round(mean_admissions, digits = 2), 
                              "<br>",
                              "2018-2019 mean admissions: ", 
                              round(mean_20182019_admissions, digits = 2)))) +
  
  geom_line(aes(x = week_ending,
                 y = mean_admissions,
                colour = simd_quintile,
                group = simd_quintile)) +
  scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") +
  scale_y_sqrt() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size =7)) +
   geom_vline(xintercept = as.numeric(as.Date("2020-01-01")), linetype=4, colour = "grey50")+
  geom_vline(xintercept = as.numeric(as.Date("2021-01-01")), linetype=4, colour = "grey50")+
  geom_vline(xintercept = as.numeric(as.Date("2022-01-01")), linetype=4, colour = "grey50")+
  labs(title = "Mean admissions per SIMD \n",
       x = "Date",
       y = "Mean admissions",
       colour = "SIMD")

`summarise()` has grouped output by 'hb_name', 'week_ending', 'simd_quintile'. You can override using the `.groups` argument.Warning: Ignoring unknown aesthetics: text

p2 %>% 
  ggplotly(tooltip = "text") %>% 
  config(displayModeBar = FALSE)